body .main-container {
  max-width: 100% !important;
  width: 100% !important;
}
body {
  max-width: 100% !important;
}

body, td {
   font-size: 16px;
}
code.r{
  font-size: 14px;
}
pre {
  font-size: 14px
}

Let us set some global options for all code chunks in this document.

knitr::opts_chunk$set(
  message = FALSE,    # Disable messages printed by R code chunks
  warning = FALSE,    # Disable warnings printed by R code chunks
  echo = TRUE,        # Show R code within code chunks in output
  include = TRUE,     # Include both R code and its results in output
  eval = TRUE,       # Evaluate R code chunks
  cache = FALSE,       # Enable caching of R code chunks for faster rendering
  fig.align = "center",
  out.width = "100%",
  retina = 2,
  error = TRUE
)
#klippy::klippy()
rm(list = ls())
set.seed(1982)

1 Preprocessing

Let us now load some required libraries.

# Load required libraries
library(INLA)
library(inlabru)
library(rSPDE)
library(MetricGraph)

library(plotly)
library(dplyr)
library(tidyr)
library(MASS)
library(glmnet)
library(car)
library(ggplot2)
library(plotly)

library(sf)

library(here)

We now define some R functions with self-explanatory names.

# Functions to manipulate the data
standardize <- function(vector) {return((vector - mean(vector)) / sd(vector))}
minmax_scaling <- function(x) {return((x - min(x)) / (max(x) - min(x)))}
convert_to_binary <- function(input_vector) {return(ifelse(input_vector != 0, 1, 0))}

1.0.1 Prepare the data

The following chunk uses the function creates_covariates_on_mesh() defined on Preprocessing/28Window.case.file4.addcovariatesonmesh.R to build a mesh (on the same graph as the one used in this file) and compute the value of all available covariates on the mesh nodes.

Go to Preprocessing 8.0.1 Raw code to see the body of the creates_covariates_on_mesh() function.

h = 0.05
source(here("Preprocessing/28Window.case.file4.addcovariatesonmesh.R"))
creates_covariates_on_mesh(h)
Sys.sleep(5)

################################################################################
################################# PREPARE THE DATA #############################
################################################################################

# loading the data
load(here("Data_files/data_on_graph_with_covariates_no_consecutive_zeros_partialtomtom.RData"))
load(here("Graph_objects/graph_construction_28_03_2024partialtomtomwhichlonglatsf.RData"))
load(here("Data_files/data_on_mesh_with_covariates_partialtomtom.RData"))


data = data_on_graph_with_covariates %>%
  mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>%
  mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
  mutate(across(c("SpeedLimit", "density_per_hour"), standardize))


mesh = data_on_mesh_with_covariates %>%
  mutate(across(starts_with(c("class_", "upto")), list(ind = convert_to_binary))) %>% # this creates new columns
  mutate(across(c("bus", "signal", "stop", "crossing"), ~round(., 5))) %>%
  mutate(across(c("SpeedLimit", "density_per_hour"), standardize))

1.0.2 Add data and build mesh

Let us add the speed observations to the graph and build the mesh. Observe that mesh already has the value of the available covariates on the mesh nodes.

sf_graph$add_observations(data = data, group = "day", clear_obs = TRUE)
sf_graph$build_mesh(h = h)
Sys.sleep(5)

2 Modeling

2.0.1 Stationary model

  • Observe that we are considering replicates.
stat.time.ini <- Sys.time()
################################################################################
################################# STATIONARY MODEL #############################
################################################################################

rspde_model_stat <- rspde.metric_graph(sf_graph,
                                         parameterization = "matern",
                                         nu = 0.5)

data_rspde_bru_stat <- graph_data_rspde(rspde_model_stat,
                                        repl = ".all",
                                        loc_name = "loc")

cmp_stat = speed ~ -1 +
  Intercept(1) +
  SpeedLimit +
  #density_per_hour +
  #bus +
  #signal +
  #stop +
  #crossing +
  #upto1_ind +
  #bus_number +
  field(loc, model = rspde_model_stat,
        replicate = data_rspde_bru_stat[["repl"]])

rspde_fit_stat <-
  bru(cmp_stat,
      data = data_rspde_bru_stat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )
stat.time.fin <- Sys.time()
print(stat.time.fin - stat.time.ini)
## Time difference of 35.79432 secs
fit.rspde = rspde.result(rspde_fit_stat, "field", rspde_model_stat)

# summaries
# summary(rspde_fit_stat)

rspde_fit_stat$summary.fixed
rbind(rspde_fit_stat$summary.hyperpar,
      rspde_fit_stat$summary.hyperpar[1,]^(-0.5),
      summary(fit.rspde))
c(unlist(rspde_fit_stat$dic[c("dic", "dic.sat", "p.eff")]),
  unlist(rspde_fit_stat$waic[c("waic", "p.eff")])) |> t() |> as.data.frame()
rspde_fit_stat$mlik[2,1]
## log marginal-likelihood (Gaussian) 
##                          -55592.19

2.0.2 Nonstationary model

  • Observe that we are using the computed parameters from the stationary model as initial values for the nonstationary models.
nonstat.time.ini <- Sys.time()
################################################################################
############################# NON STATIONARY MODEL #############################
################################################################################

B.sigma = cbind(0, 1, 0, mesh$SpeedLimit, 0)
B.range = cbind(0, 0, 1, 0, mesh$SpeedLimit)

rspde_model_nonstat <- rspde.metric_graph(sf_graph,
                                            start.theta = c(fit.rspde$summary.log.std.dev$mode,
                                                            fit.rspde$summary.log.range$mode,
                                                            rep(0, 2)),
                                            B.sigma = B.sigma,
                                            B.range = B.range,
                                            parameterization = "matern",
                                            nu = 0.5)

data_rspde_bru_nonstat <- graph_data_rspde(rspde_model_nonstat,
                                        repl = ".all",
                                        loc_name = "loc")

cmp_nonstat = speed ~ -1 +
  Intercept(1) +
  SpeedLimit +
  #density_per_hour +
  #bus +
  #signal +
  #stop +
  #crossing +
  #upto1_ind +
  #bus_number +
  field(loc, model = rspde_model_nonstat,
        replicate = data_rspde_bru_nonstat[["repl"]])

rspde_fit_nonstat <-
  bru(cmp_nonstat,
      data = data_rspde_bru_nonstat[["data"]],
      family = "gaussian",
      options = list(verbose = FALSE)
  )

nonstat.time.fin <- Sys.time()
print(nonstat.time.fin - nonstat.time.ini)
## Time difference of 37.44091 secs
# summaries
# summary(rspde_fit_nonstat)

rspde_fit_nonstat$summary.fixed
rbind(rspde_fit_nonstat$summary.hyperpar,
      rspde_fit_nonstat$summary.hyperpar[1,]^(-0.5),
      summary(rspde.result(rspde_fit_nonstat, "field", rspde_model_nonstat)))
c(unlist(rspde_fit_nonstat$dic[c("dic", "dic.sat", "p.eff")]),
  unlist(rspde_fit_nonstat$waic[c("waic", "p.eff")])) |> t() |> as.data.frame()
rspde_fit_nonstat$mlik[2,1]
## log marginal-likelihood (Gaussian) 
##                           -55598.8
save.image(here(paste0("Models_output/", rmarkdown::metadata$title, ".RData")))